|
The ziggurat algorithm is an algorithm for pseudo-random number sampling. Belonging to the class of rejection sampling algorithms, it relies on an underlying source of uniformly-distributed random numbers, typically from a pseudo-random number generator, as well as precomputed tables. The algorithm is used to generate values from a monotone decreasing probability distribution. It can also be applied to symmetric unimodal distributions, such as the normal distribution, by choosing a value from one half of the distribution and then randomly choosing which half the value is considered to have been drawn from. It was developed by George Marsaglia and others in the 1960s. A typical value produced by the algorithm only requires the generation of one random floating-point value and one random table index, followed by one table lookup, one multiply operation and one comparison. Sometimes (2.5% of the time, in the case of a normal or exponential distribution when using typical table sizes) more computations are required. Nevertheless, the algorithm is computationally much faster than the two most commonly used methods of generating normally distributed random numbers, the Marsaglia polar method and the Box–Muller transform, which require at least one logarithm and one square root calculation for each pair of generated values. However, since the ziggurat algorithm is more complex to implement it is best used when large quantities of random numbers are required. The term ''ziggurat algorithm'' dates from Marsaglia's paper with Wai Wan Tsang in 2000; it is so named because it is conceptually based on covering the probability distribution with rectangular segments stacked in decreasing order of size, resulting in a figure that resembles a ziggurat. ==Theory of operation== The ziggurat algorithm is a rejection sampling algorithm; it randomly generates a point in a distribution slightly larger than the desired distribution, then tests whether the generated point is inside the desired distribution. If not, it tries again. Given a random point underneath a probability density curve, its ''x'' coordinate is a random number with the desired distribution. The distribution the ziggurat algorithm chooses from is made up of ''n'' equal-area regions; ''n'' − 1 rectangles that cover the bulk of the desired distribution, on top of a non-rectangular base that includes the tail of the distribution. Given a monotone decreasing probability density function ''f''(''x''), defined for all ''x'' ≥ 0, the base of the ziggurat is defined as all points inside the distribution and below ''y''1 = ''f''(''x''1). This consists of a rectangular region from (0, 0) to (''x''1, ''y''1), and the (typically infinite) tail of the distribution, where ''x'' > ''x''1 (and ''y'' < ''y''1). This layer (call it layer 0) has area ''A''. On top of this, add a rectangular layer of width ''x''1 and height ''A''/''x''1, so it also has area ''A''. The top of this layer is at height ''y''2 = ''y''1 + ''A''/''x''1, and intersects the density function at a point (''x''2, ''y''2), where ''y''2 = ''f''(''x''2). This layer includes every point in the density function between ''y''1 and ''y''2, but (unlike the base layer) also includes points such as (''x''1, ''y''2) which are not in the desired distribution. Further layers are then stacked on top. To use a precomputed table of size ''n'' (''n'' = 256 is typical), one chooses ''x''1 such that ''x''''n'' = 0, meaning that the top box, layer ''n'' − 1, reaches the distribution's peak at (0, ''f''(0)) exactly. Layer ''i'' extends vertically from ''y''''i'' to ''y''''i''+1, and can be divided into two regions horizontally: the (generally larger) portion from 0 to ''x''''i''+1 which is entirely contained within the desired distribution, and the (small) portion from ''x''''i''+1 to ''x''''i'', which is only partially contained. Ignoring for a moment the problem of layer 0, and given uniform random variables ''U''0 and ''U''1 ∈ [0,1), the ziggurat algorithm can be described as: # Choose a random layer 0 ≤ ''i'' < ''n''. # Let ''x'' = ''U''0''x''''i''. # If ''x'' < ''x''''i''+1, return ''x''. # Let ''y'' = ''y''''i'' + ''U''1(''y''''i''+1 − ''y''''i''). # Compute ''f''(''x''). If ''y'' < ''f''(''x''), return ''x''. # Otherwise, choose new random numbers and go back to step 1. Step 1 amounts to choosing a low-resolution ''y'' coordinate. Step 3 tests if the ''x'' coordinate is clearly within the desired density function without knowing more about the y coordinate. If it is not, step 4 chooses a high-resolution y coordinate, and step 5 does the rejection test. With closely spaced layers, the algorithm terminates at step 3 a very large fraction of the time. Note that for the top layer ''n'' − 1, however, this test always fails, because ''x''''n'' = 0. Layer 0 can also be divided into a central region and an edge, but the edge is an infinite tail. To use the same algorithm to check if the point is in the central region, generate a fictitious ''x''0 = ''A''/''y''''1''. This will generate points with ''x'' < ''x''''1'' with the correct frequency, and in the rare case that layer 0 is selected and ''x'' ≥ ''x''''1'', use a special fallback algorithm to select a point at random from the tail. Because the fallback algorithm is used less than one time in a thousand, speed is not essential. Thus, the full ziggurat algorithm for one-sided distributions is: # Choose a random layer 0 ≤ ''i'' < ''n''. # Let ''x'' = ''U''0''x''''i'' # If ''x'' < ''x''''i''+1, return ''x''. # If ''i'' = 0, generate a point from the tail using the fallback algorithm. # Let ''y'' = ''y''''i'' + ''U''1(''y''''i''+1 − ''y''''i''). # Compute ''f''(''x''). If ''y'' < ''f''(''x''), return ''x''. # Otherwise, choose new random numbers and go back to step 1. For a two-sided distribution, of course, the result must be negated 50% of the time. This can often be done conveniently by choosing ''U''0 ∈ (−1,1) and, in step 3, testing if |''x''| < ''x''''i''+1. 抄文引用元・出典: フリー百科事典『 ウィキペディア（Wikipedia）』 ■ウィキペディアで「ziggurat algorithm」の詳細全文を読む スポンサード リンク
|